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Abstract The aim of this work is to consider multiscale algorithms for solv- 
ing PDEs with Galerkin methods on bounded domains. We provide results 
on convergence and condition numbers. We show how to handle PDEs with 
Dirichlct boundary conditions. We also investigate convergence in terms of the 
mesh norms and the angles between subspaces to better understand the differ- 
ences between the algorithms and the observed results. We also consider the 
issue of the supports of the RBFs overlapping the boundary in our stability 
analysis, which has not been considered in the literature, to the best of our 
knowledge. 



1 Introduction 

Radial basis functions (RBFs) have been increasingly important in the area 
of approximation theory. For solving partial differential equations (PDEs) 



RBFs have been studied more for meshlcss collocation ( Gicsl fc Wendland 



20071 lFasshauerl2007l) however their use with Galerkin methods have also been 



considered ( Wendland 1998 fclfc ). Two excel lent r ecent book s cover ing practical 



and theoretical issues are Fasshauer ( 2007 ) and Wendlandl (|2005l ). A function 
<P : W l -+ M. is said to be radial if there exists a function <f> : [0, oo) — > ffi such 
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that <f>(x) = 0(||x|| 2 ) for all x G R d , where 



denotes the usual Euclidean 



norm in R . Then we can define an RBF for a given centre Xj € R as 

#(xi) :=0(||x-Xi|| 2 ). 

A practical issue that arises is that of which scale to use for the radial basis 
functions. A small scale will lead to a sparse and consequently well-conditioned 
linear system, but at the price of poor approximation power. Conversely, a 
large scale will have better approximation power but at the price of an ill- 
conditioned linear system. 

Many examples may naturally exhibit multiple scales, for example, con- 
structing an approximation for the height of the earth's surface may suggest 
a "large scale" to be used over desert regions and a "fine scale" over areas of 
high variability, such as the Himalayas. It appears much more appropriate to 
allow different scales in different regions. Of course, this comes at the price of 
having to select which scales to use in which regions but this is not the topic 
of this paper. 

The multiscale algorithms proposed in this paper are constructed over mul- 
tiple levels, in which the residual of the current stage is the target function for 
the next stage, and at each stage, RBFs with smaller support and with more 
closely spaced centres will be used as basis functions. 

Such a multiscale algorit hm for interpolation was first proposed in Floater fc Iskel 



(|1996l ) and ISchabackl (|1996f ) but without any theoretical grounding. Theoret- 
ical con vergence was prove n in the case of the data points being located on a 
sphere (|Le Gia et alJl2010h a nd then extende d to interpolation and approxi- 
matio n on bounded dom ains ( Wendlandl 2010h . 

In I Wendlandl ()1998a) we can find various experiments with two multiscale 
algorithms for constructing Galerkin approximations to PDEs on bounded 
domains. These two algorithms are studied in this paper. There was no proof of 
convergence for the first algorithm which we call the multiscale algorithm. For 
the second multiscale algorithm, which we call the nested multiscale algorithm, 
convergence was proven, making use of the fact that the weak formulation of 
a PDE can be interpreted as a Hilbert space projection method. Numerical 
experiments with both multiscale algorithms were also given. 

The aim of this work is to consider multiscale algorithms for solving PDEs 
with Galerkin methods on bounded domains. We provide results on conver- 
gence and condition numbers. We show how to handle PDEs with Dirichlet 
boundary conditions which was not covered in the earlier papers. We also 
investigate convergence in terms of the mesh norms and the angles between 
subspaces to better understand the differences between the two algorithms 
and the observed results. We also consider the issue of the supports of the 
RBFs overlapping the boundary in our stability analysis, which has not been 
considered in the literature, to the best of our knowledge. 

Related multileve l methods for the numerical solution of PDEs include do- 
main decomposition ( Smith et al.lll996h and the use of smoothing as a postcon- 
ditioner with an iterative ( multilevel) method (lFasshauerlll999l) . The hetero- 
geneous multiscale method (jWeinan et al.ll2007l ) differs in that it is multiscale 
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in the time domain, whilst we work with m ultiple scales in the spatial domain. 
Multigrid methods ( Brenner fc Scott 20081) require the use of a computational 
mesh, whilst our approach is meshfree. 

In the next section we review the background material required for the 
multiscale algorithms. Sections [3] and [4] describe how to construct Galerkin 
solutions to PDEs using radial basis functions for PDEs without and with 
Dirichlct boundary conditions respectively. Sections [5] and |5] consider the mul- 
tiscale and nested multiscalc algorithms. Section [7] provides the results of sev- 
eral examples using the two multiscale algorithms. Finally Section [5] provides 
an analysis of the convergence of the two algorithms. 



2 Preliminaries 

In this paper, we will use (scaled) compactly supported radial basis functions 
to construct multiscale approximate solutions to PDEs, that is, we form the 
solution over multiple levels. We will work with a given domain fl C M. d . A 
kernel <P : fl x fl — > K is also given. 

At each level, we will have a finite point set X C fl. We will define the fill 
distance as 

hx,n ■= sup min ||x - Xj|| 2 , 

which is a measure of the uniformity of the points in X with respect to fl. Then 
for example, at each level i, we denote the fill distance by hi- The selection of 
point sets with fill distances decreasing in a specific way will form one of the 
requirements for convergence of our algorithms. 

The functions that we will be concerned with are defined on a bounded 
domain fl with a Lipschitz boundary. As a result, there is an extension operator 
for func tions defined in Sobol ev spaces which is presented in the following 
lemma ( Brenner fc Scott 120081 Theorem 1.4.5). For further details, we refer 
the reader to lSteinl (|l970l) . 



Lemma 2.1 Suppose fl C R d has a Lipschitz boundary. Then there is an 
extension mapping E : H T (fl) — > H T (R d ), defined for all non-negative integers 
r, satisfying Ev\q = v for all v € H T (fl) and 

\\Ev\\ H r {Rd) < C\\v\\ H , {n) . 

In this paper, C will denote a generic constant. 

Since we also have || f|| ( r?) < ||^ w ll_ff T (R ci )i this means that when we need 
to consider the H T (fl) norms of the errors at each level, we can carry out our 
error analysis in the ff r (M d )-norm. This is advantageous, since we then have 
for g E H T (R d ) 



\ 2 H* 



|^)| 2 (l + ||u;||2) r du;, (1) 
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upon defining the Fourier transform as 

= (2^)~ d/2 / g(x)e-" 



At each level, we will also require a scaled version of the kernel <P : Q x J? — >• M. 
For our unseal ed kernel we wil l use a Wendland compactly supported radial 
basis function (|Wendlandll2005l) . With a (level-specific) scaling parameter 5 > 
0, we can define the scaled kernels as 



*«(x,y) :=5~U 



y a 



(2) 



Appropriate selection of the scaling parameters will also prove to be one of 
the important ingredients for convergence of our multiscale algorithms. 

For the Wendland basis functions, the re exist two cons tants < c\ < c 2 
such that their Fourier transforms satisfy (|Wendlandll2005l) 



a 1 



< < c 2 (1 + \\u\\l) 



<jJ G 



(3) 



where r = (d + 2k + l)/2 where d is the spatial dimension an d fc is the 
smoot hness parameter. For further details, we refer the reader to Wendlandl 
(|2005t ). The native space Af<p(R d ) of # consists of all functions g G L 2 {M. d ) such 
that 



-du> < oo. 



(4) 



Taking ([T]) and ^ together shows that A/"<z>(R d ) is norm-equivalent to the 
Sobolev space H T (R d ). 

Consequently the Fourier transform of <P$, ^(w) = <P(8u)), satisfies 



ci (l + 5 2 \\u\ 



< < c 2 (l + <J 2 H||) 



W G 



(5) 



We will need norm equivalence as stated in the following lemma. 



Lemma 2.2 For every 6 <G (0, S a ] and for all g G H T (R. d ), there exist con- 
stants < C3 < C4 such that 

C3\\gh s < h\\Hum") < c4<5" T ||ff||<f,- 



Proof The case S a < 1 was proven in Wendlandl ( 2010h with C3 = c^ 2 and 



Ci = 

have 



.1/2 



To extend this to the case where S a > 1, note that for 5 > 1 we 



(1 + |M| 2 ) T = s-^ (5 2 + * 2 |M| 2 ) T > (1 + §2 ]Mt) 
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Together with (JXJ) , ([5]) and Q we can see that 

IMIff T (R d ) 



- ' ^(u>)\ 2 (1 + \\u:\\lY du: 



>S- 2T J |5M 2 (l + 5 2 |M|^du; 



2r f \9("V 2 



MPa + iMjsydu; 



>ci8-^\\g\\% s . 

For the lower bound, we can just use 6 > 1 directly to derive 

< J \g(u;)\ 2 (l + S 2 \\uj\\l) T duj 

< C2\\g\\l 5 , 

using ((5|) and (HJ). Then setting C3 := c^ 2 min(l, 6~ T ) and c 4 := c^ 2 max(l, SJJ 
completes the proof. □ 



3 PDEs with Neumann and/or Robin boundary conditions 

In this section, we consider a second order PDE which has homogeneous Neu- 
mann and/or Robin boundary conditions. For example, such a PDE with Neu- 
mann boundary conditions is given by 

Cu = f in i?, (6a) 

— = on dSi, (6b) 
on 

where £ is a second order elliptic differential operator, n denotes the outward 
unit normal vector and dfi denotes the boundary. The weak formulation is 
given by 

a(u, v) = (/, v) L2(n) V v G V, (7) 

where V = H 1 (f2). We assume that C and / are such that a(u,v) is a strictly 
coercive and continuous bilinear form defined on V X V and (/, w)i 2 (r2) ^ s a 
continuous linear form defined on V . By the Lax-Milgram theorem, ([7]) has a 
unique solution u gV. We will also require u G H 2 {fi) with spatial dimension 
d < 3. 
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Galcrkin approximation seeks to solve Q with a finite dimensional sub- 
space Vn C V. In other words, the Galerkin approximation un is the solution 
of 

u N £ Vn ■ a(u N ,v) = (f,v) L ^ n) Vu G Vn- (8) 

We will consider fl to be a bounded domain with a Lipschitz boundary, which 
means that we can apply the extension operator to use norms in K d as ex- 
plained in the introduction. For further information on we ak formulation of 
PDEs and Galerkin approximation, we refer the reader to iBrenner fc Scott 
(|2008h . 

Since the PDE does not have Dirichlet boundary conditions, we can use 
the entire Sobolev space H 1 {Q) rather than the subspace £T 1 (f2) consisting 
of functions with zero boundary values which can occur with pure Dirichlet 
boundary conditions. 

We will consider finite dimensional subspaces Vn Q V of the form 

V N :=span{#(-, Xj ) : 1 <j< N} 

where <S> : M. d — > R is at least a C 1 -function and there are N centres { Xj : 1 < 
j < N}. In this case our approximation takes the form 

N 

u N =^]c^(-,x J ), 

i=i 

and the weak formulation with this approximation given by 

a{u Nl v) = (/, v) L2(n) y v £ Vn, 
results in a linear system 

Ac = f 

where the entries of the stiffness matrix are given by 

Atj = o(#(-,x j ),*(. > 3c,-)) ( 9 ) 

and 



fi = J /(x)<P(x,x 4 )dx. 



We have the following result from IWendlandl (|l998Zft . 

Theorem 3.1 If u £ H 2 (Q),d < 3, is the solution to the variational problem 
and un £ Vn is the solution of ([8|. where Vn is generated with X sat- 
isfying h < ho for hg small enough and fixed, then the error can be bounded 
by 

unWhHs-2) <Ch\\u\\ H 2 (sl) . 
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Lemma 3.2 Consider the adjoint variational problem 

a(v,u) = (f,v) LAn) VveV. (10) 
If the solution u satisfies the regularity estimate 

\\u\\ HH n) < C\\f\\ L2m (11) 

for both the variational problem ([7]) and the adjoint variational problem (|10[) . 
then we have the following error bound 

\\u - wjv||i 2 (fi) < Ch\\u - u N \\ H i {n) < Ch 2 \\u\\ H 2 (n) . 

Proof Wc follow a duality argument as in ([Brenner k, Scott 2008. Theorem 
5.7.6). With un as defined in (|8]), let w be the solution to the adjoint problem 

a(v,w) = (u - u N ,v) L2{n) VveV, 

and let the Galerkin approximation be given by w. Then since the bilinear 
form a(-, •) is bounded, u £ H 2 (f2), and with Theorem 13.11 we have 

ll u -"Jv||i 2 ( fi) = (u-u N ,u-u N ) L2(n) 

= a{u — un, w) = a(u — un, w — w) 

< C\\u - Ujv||j?i(fi) \\w - w\\ H i(n) 

< Ch\\u - ujvllffi(fi) \\w\\h^-(S2) 

< Ch\\u - u N \\ H i (n) \\u - UnWl^q), 

where in the second last line, we use the regularity estimate (jllj) and the result 
follows with another application of Theorem 13.11 □ 



We note that is known to hold (jBrenner fc Scottll2008l p. 139) 

— if f2 has a smooth boundary and the problem has pure Dirichlet or pure 
Neumann boundary conditions; 

— if d = 2 and J? is convex and the problem has pure Dirichlet or pure 
Neumann boundary conditions. 



4 PDEs with Dirichlet boundary conditions 

In this section we will consider a PDE with Dirichlet boundary conditions and 
we seek error estimates similar to those in the previous section. For ease of 
description, wc will consider the following boundary value problem 

-Au = f in (12a) 
u = g on dn, (12b) 



wh ere A denotes the L aplacian in R d . 

Nitschel (|l970-197lh proposed minimising the functional J(v — u) where 
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JO) := / |Vtor - 2 / w (Vro ■ n) + p N J w, 
n an an 

for all v £ Vn with u being the solution of (|12p and where n denotes the 
outward unit normal vector and V the gradient operator. The parameter (3n > 
depends only on the subspace Vn- The approximation un is given by J{un — 
u) := infugy^ J(v — u). With / and g from (fP2|) . we can compute Un since 

J{v -u) = J(v) + J{u) - 2ljfv + jg {fi N v - Vv ■ n) 

\n an 



Then the variational form to approximate ([12)) becomes: find un € Vat such 
that for all v £ Vn 

a D {u N ,v) = i D {v), (13) 
where we use the subscript D to denote the Dirichlet boundary conditions and 

cld(u, v) := Jvu-Vv — J v (Vu • n) — J u(yv-n)+(3N J uv (14a) 
n an an an 

io{v) := / fv- J g (Vv -n) + p N J vg. (14b) 

n an an 

It can be shown that the variational form using Nitsche's method leads to 
variational consiste ncy, in the sense that if u is sufficiently regular, then 



(jBlowev et al.ll2003l p. 119) 

a D (u,v) = Id{v), V v £ V N - 
If there exists a positive constant Cn such that 

||Vw • n|| L2(ar2) < ^L\\Vv\\ L2(n) Vv £ V N , (15) 

with 8 being the support of the radial basis functions, then selecting 

Pn = j (16) 

with C5 > 2C N will ensure that the bilinear form cld{-, •) is symmetric positive 
definite. This choice of C5 will also ensure that an is coercive since 

a D (v,v) = \\Vv\\l 2[n) - 2 J v {Vv ■ n) + P N \\v\\l 2idQ) 

an 

> \\^ v \\l 2 (n) ^ jq^M\L 2 (an)\\^v\\L 2 (n) + Pn\\v\\ 2 L2(9[2) 

>\\m\i 2{Q) + (^ N - 2 -^) \\v\\i 2(an) 

>C\\vf m{n) ,\/v £ Vn 
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where we have used the Cauchy-Schwarz and Friedrichs inequal ities. (| 1 5J) and 
that 2xy < x 2 + y 2 . We recall that the Friedrichs inequality (|Maz'vall201ll ) 
states that 

' u\ 2 <C 




for fl being a bounded domain for which the Gauss-Green formula holds. 
Continuity follows since 

\a D (u,v)\ < \u\ H i {n) \v\ H i( n) + C N /5 (\\v\\ L2{dn) \\\7u\\ L2{n) + 

\\u\\L 2 (dO)\\^v\\ L2 ^n)) + PN\\u\\L 2 (dn)\\v\\L 2 (dn) 

< \u\ H i(n)\v\ H i {n) + C (\\v\\ L . 2{n) \\Vu\\ L2{n) + 
lklU 2 (fi)||Vv||z, 2 (fi)) + PN\\u\\ L2{n) \\v\\ L2(n) 

< C\\ u \\m(i:2)\\v\\m({2), Vu,u e V N 

where we have used the Cauchy-Schwarz inequality, (|15[) and the Sobolcv trace 
embedding theorem. 

Nitsche also proved that the optimal error estimates of Theorem 13.11 and 
Lcmma [3 . 2l hold in this setting if, in addition to the requirement of selecting /3/v 
satisfying ([To]) , there exists a s u £ Vjv such that for u £ H 2 (f2), the following 
error bounds hold for k G {0, 1} 

\\u - s u \\ Hk(n) < Ch 2 ~ k \\u\\ H 2 {n) , (17a) 
\\u - s u \\ HHen) < Ch 3 / 2 - k \\u\\ H2(n) . (17b) 

With our choice of co mpactly supported radial basis functions, this require- 
ment is known to hold (IWendlandlll998^ . 



Note that the most challenging aspect of Nitsche's method is the derivation 
of the weak form and the selection of the stabilisation parameter (3jy. Both the 
weak form and the choice of the parameter ^jv depend on the PDE as well as 
the Dirichlet boundary conditions. 



5 Multiscale Galerkin approximation 

In this section, we consider a multiscale algorithm for constructing a Galerkin 
approximation where we use the residual from the previous level as the target 
for each subsequent level. We define the approximation at level i as u,; := 
ujV; with centres iVj and the approximation space at level i as Vi := V/\r. 
The algorithm is given in Algorithm [T] The bilinear form a(-, •) used in this 
algorithm is the unmodified bilinear form in the case of a PDE with Neumann 
or Robin boundary conditions and the Nitsche's method bilinear form ao(-, ■) 
in the case of a PDE with Dirichlet boundary conditions. 

The algorithm as stated uses the same bilinear form at each level and it 
is the approximation space Vi which changes. However the Nitsche's method 
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Algorithm 1: Multiscale Galerkin approximation 
Data: n: number of levels 

{Xi}™ =1 : the set of nested centres for each level i, with mesh 
norms at each level given by hi satisfying cfihi < /ij+i < fxhi 
with fixed /i G (0, 1), c G (0, 1] and h\ sufficiently small 
{3i}i=i '■ * ne sca le parameters to use at each level, 
satisfying 5{ = vhi, v a fixed constant. 

begin 

Set uq = 

for i = 1,2, ... ,n do 

With the level-specific approximation subspace 

Vi := span{#5 4 (•, x), x G X,} solve the Galerkin approximation 

given by 

Find Si G Vi : a(sj,u) = (f,v) L2 ^ 2) - a(Ui-i,v) \fv G 
Update the solution according to 

= + Si 



Result: Approximate solution at level n, u. 
The error at level n, e„ :— u — u n . 



bilinear form Od(-, •) will vary at each level since the value of is proportional 
to S^ 1 and hence so is /3/v- This means that we will need to select the value of 
/3jv corresponding to the last level and to use this for all previous levels. This 
will also mean that we will need to know the number of levels in advance. 

Henceforth we will simply refer to the bilinear form as a(-, ■). This should 
cause no confusion as we have the same error bounds in both cases, as well 
as coercivity and continuity, and the multiscale algorithm follows the same 
steps in both cases. We require one more lemma before we can analyse the 
convergence of the multiscale algorithm. 

Lemma 5.1 Let fl C R d be a bounded domain with a Lipschitz boundary. 
Then for Algorithm]^ and for a given level i > 1, we have the following bound 
on the H l (Q) error between subsequent levels. 



||ei||jyi(n) < C'||ej_x||jii(fl), 



where e, is defined in Algorithm]^ 
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Proof We will firstly show that Sj is the Galerkin approximation of e,_i. We 
have 

a(si,w) = {f,w) L2{n) - a(ui-i,w), w e V t 
= a(u, w) — a(v,i-i,w) 
= a(u — Ui-i, w) 
= a(ei-i,w), 

where we have used the variational form of the PDE and the linearity in the 
first argument of the bilinear form a(-, •). Hence on setting w = s, we obtain 

a(e,_i - Si) = 0. 

Upon noting that e; = e^—i — Sj, it follows easily that 

a(e l , e») = a(e,_i, e ?; _i) - a(sj, Sj), 

and since the bilinear form a is continuous and coercive 

IKIIfl- 1 ^) + Ikillff 1 ^) - C|| e i-1 II ffi(fT)! 

from which the result follows. □ 

The following theorem and corollaries are our main results for the conver- 
gence of the multiscale Galerkin approximation. For the error analysis, we will 
need the norm 




F:-< 

As in Lemma this norm satisfies 

C7\\u\\^ j < \\u\\ H i {Rd) < cs^lHlsy (18) 

Theorem 5.2 Let fl C M. d be a bounded domain with a Lipschitz boundary. 
Then for Algorithm]^ there exists a constant ot\ > such that 

||-Eej||tfi + i < Qiill^ej-ill^ for i = l,2,... 

where Eei is the extension operator defined in Lemma \2.1\ applied to e; . The 
constant a>i satisfies ct\ < 1 if in Algorithm\^v is sufficiently small and (i is 
sufficiently large. 

Proof Using ([4]) and ([5]), we can write 

||^|||. +1 < - ( |^M| 2 (1 + <5? +1 |M|2) du; =: - (h + I 2 ) 

Cl J Ci 

R d 

with 

ll-l|2<^ 
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h~ J iBeiHI^l + ^HlDdw. 



Now wc consider the first integral where we can use that <5j_|_i||u>||2 < 1 and 
then Lemmas 13.21 and (1181) . 



h < 2 J |-Ee;(w)| 2 &u) 



\\u\h<Tr— 



< 2 ||^||^ (Rd) < C |N|| 2( , 2) < C\ 2 Heill^^) 

< C/i? Hei-iH^^) <c(|) ll^-ill* 
= C^- 2 || J Be i _ 1 ||| i , 

where we have also used Lemma l5.ll For J2, since <5j_|_i||ci>||2 > 1, we have that 

(1 + <5. 2 +1 |M| 2 ) < 2<5 2 +1 ||u,|| 2 < 2<5 2 +1 (l + ||u, 2 ) . 
Then again using Lemma 15.11 shows that 



h < 25 2 +1 ||Se i || 2 i - 1(Rd) < C8l +1 \\ei\\ 2 H i (n) 



<C8i + Aei-i\\l, {n) <c{^j WEe^g 
< CfWEei-iWl.. 



(19) 



Combining our results for I\ and 1% and now writing C\ and C2 for the two 
constants appearing in the bounds of the expressions for I\ and I2 respectively, 
we have that 

||£ei||! i+1 < (v- 2 Cx/ Cl + /x 2 C 2 / Cl ) ||£ei_i||| 4: 
and the result follows with 

ai := (i^CiM + M 2 C 2 / Cl ) 1/2 . (20) 

□ 

Corollary 5.3 There exists a constant C3 > suc/i that 

\\u-u n \\ L2 ^ 2) <C 3 a1\\u\\ H i^ 2 ) for n=l,2,... (21) 

T/ims u n resulting from Algorithm]]] converges linearly to u in the L2-norm if 
ai < 1. 
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Proof Using Lemmas 13.21 and 12.21 again, we can see that 

\\u-u n \\ L ^n) = ||e„||z, 2 (fl) < Ch n \\e n \\ H i(a) 

< Ch n \\Ee n \\ H i(ui) < Ch n 5~\ x \\Ee n \\ 9n+1 

< C\\Ee n \\v n+l , 

since 

far). hn. - 1 



5 n+ i vh n+1 cfiv 
Now we can apply Theorem 15.21 n times, and noting that uq = 0, leads to 

\\u-u n \\ L2(i2) < Ca1\\Eu\\ Wl < Ca1\\Eu\\ m{Rd) < Ca^\\u\\ H i {n) . 

□ 

5.1 Condition number 

In this subsection, we present upper and lower bounds for the eigenvalues of 
the multiscale Galerkin algorithm. Since the Galerkin approximation matrix is 
symmetric and positive definite, we know that the condition number is given 

by 

-(A) = (22) 

where A max (A) and A m i n (A) denote the maximum and minimum eigenvalues 
of A as defined in ©. 

Theorem 5.4 Let <P be a positive definite kernel generating H T (R d ) with 
t > d/2. Let <S>i := <?(-,x,;) and assume that there exists a constant c§ > 
depending only on <P and Q such that 

a T (F-c 6 G)a>0, V a e R w , (23) 

which means that F — cgG is positive semi-definite and where 

Fi,j = (*i.*j> ff i(fl). 
Then the condition number of A can be bounded by 

, -. \ 4-r-2d 

k(A) < C - 
where the constant C is independent of the point set X. 
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Proof Since v = ^2 i=1 oti&i and with with the coercivity of a(-, •) and (|2"3"]l . wc 
have that 

a(v,v) >C\\v\\ 2 HHn) =Ca T Fa T 

> Ca T Ga > CX min (G)\\a\\l 

where A m ; n (G) is the minimum eigenvalue of G. From IWendlandl (<1998dh we 
know that (@i,4>j)gimd\ is a radial function given by 



T (xi, Xj) := - (Z\0) * 0(x, - x,) + 0(x, - x 



where Z\ again denotes the Laplace op erator and * den otes convolution defined 
as f*g(x) := J f(y)g(x—y) dy. From ( Wendlandll2005 , Theorem 12.3) we know 
that we can use T to derive a lower bound on the minimum eigenvalue of G. 
Then we have that 

T(z) = (l + |z| 2 )£ 2 (z). 

From (IWendlandll2005l Theorem 10 35), we know th at $(z) > C\\z\\~ 2t and 
hence f (z) > C||z||- 4t . Then using (|Wendlandll2~b"oi Theorem 12.3) we reach 

A mi „(G) > Cq%- d . 

With the continuity of a(-, •), I Wendlandl 1998d . Theorem 14.2) and the non- 
negativity of norms, we also have the following bound on the maximum eigen- 
value 

a(v,v) < C\\v\\ 2 HHn) < C\\v\\ 2 Hl{Md) = Ca T Ga 
< CNq x d \\a\\l. 

These two bounds, in conjunction with (|2"2"]l . complete the proof. □ 

$s(-,Xi), where <P is 



We will consider (|23p further with the scaled RBFs #j : 
the C 6 Wendland function given by 

<P(x) = (l-||x||) 8 + (32||x!| 3 + 25||x|| 2 + 8| 

which is positive definite on R 2 (jWendland 2005 ). Now since the support of 
the radial basis functions is fixed, \\$i\\ h 1 ^) is fixed and is independent of 
the point set X and f2 and we can express this as 



1) 



(24) 



2tt S 




drdO. 



(25) 



o o 



In the case of the unit square, which will be used for the numerical experi- 
ments in Section [71 note that ||^j|| is minimised when the RBF centre is 
located on a corner. This can be easily verified since moving the centre in any 
direction (in ft) will keep the original area inside fl and lead to additional area 
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being inside fl and the integrand is non-negative. Then we can then bound 

Nl/fi^) b y 




(26) 

where we have also used equations from ( Fasshauer 20071 Appendix D). 
This means that 

._ Il^llg 1 (« d ) 9 

As a result, if 6 < qx, which means that there is no overlap between the vari- 
ous RBFs and F and G are diagonal, or if F and G arc diagonally dominant, 
then ([23|l will hold since then we know that positive diagonal entries will en- 
sure at least positive scmi-dcfinitcncss. Whilst we have been unable to prove 
that (|23[) holds in full generality for the unit square, it is supported by exten- 
sive numerical testing. The numerical experiments in Section [7] also provide 
empirical evidence since Corollarv l5.6l holds, which depends on Theorem 15.41 
We have the following theorem on the condition number of the multiscalc 
algorithm. 

Theorem 5.5 Let <P be a positive definite kernel generating H T (M. d ). Then 
the condition number of the Galerkin approximation matrices from Algorithm 
[7] can be bounded by 

K (A) < C , (27) 

with a constant C > independent of X and of the scaling parameter S. 

Proof At each level, we now introduce the dataset X/S = {xx/5, . . ,,x.m/S}, 
which obviously has separation radius 

qx/s - — 

and since a(-, •) is bilinear, the Galerkin approximation matrix at each level is 
A X ,5 = K^O.xO.^-.Xj))) 

rt(*(^), # (_ii)))=^ AAW . 

Then the result follows with Theorem 15.41 □ 



Corollary 5.6 If the datasets are quasi-uniform, which means that hj/qj is 
bounded above by a constant, then the condition numbers of the Galerkin ap- 
proximation matrices from Algorithm^ are bounded above by a constant. 
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Proof Algorithm[T]takcs Sj = vhj with a constant v > 1. With the assumption 
of quasi-uniformity, hj < cqj and the result follows with Theorem l5.5l □ 

We note that since we only require the bilinear form a(-, •) to be continuous 
and coercive, these theorems on the condition number will apply for PDEs with 
Robin and/or Neumann boundary conditions as well as PDEs with Dirichlet 
boundary conditions. 



6 A nested multiscale algorithm 



In this secti on we will consider another multiscale Galerkin algorithm that was 
proposed in lWendlandl (|1998a ) . This essentially extends Algorithm [T] and hence 
we can consider a PDE cither with or without Dirichlet boundary conditions. 
We refer to this as a nested multiscale algorithm because it contains inner 
and outer iterations. We will also see that this has a connection to multigrid 
methods from the finite elements literature. The details are given in Algorithm 

12 



Algorithm 2: Nested multiscale Galerkin approximation 
Data: K: number of outer levels 

n: number of inner levels 
{-X^}" =1 : the set of centres for each inner level £, with mesh 
norms at each inner level given by he satisfying 
cfiht < hi+i < nhe with fixed \i G (0, 1), c G (0, 1] and 
hi sufficiently small 

{Se}" =1 : the scale parameters to use at each inner level, 
satisfying Se = fhe, v a fixed constant. 

begin 

Set uo = 

for k = 0,1,2, K do 
for £ = 1, 2, . . . , n do 

With the level-specific approximation subspace 
Ve := span {^(-, x), x G Xe] solve for the Galerkin 
approximation given by 

Skn+e € V e : a(s kn +e,v) = {f, v ) L2 (n) ~ a (ukn+e-i, v) 
for all v G Ve. Update the solution according to 

Ukn+e = Ukn+e-1 + Skn+e 



Result: Approximate solution at level n, u n (K+i) 

The error at level n(K + 1), e n rjc+i) '■= u — u n (K+i)- 
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From Wendlandl (1998c) we have the following theorem regarding conver- 
gence. 

Theorem 6.1 Let u* denote the best approximation to u from V\ + . . . + V n 
with respect to the norm \\ ■ Wh 1 ^)- Then there exists 9 £ (0, 1) such that 

- un\\m(n) < 9 K |M|ffi(r?), 
where uk is the approximation from Algorithm^ 

Note however that this docs not mean that we have linear convergence of 
the approximation from Algorithm [5] to the true solution u. The convergence 
of the approximation from Algorithm [2] to the true solution u is given in the 
following theorem. 

Theorem 6.2 Let Q C W 1 be a bounded domain with Lipschitz boundary. Let 
<P be a kernel generating H 1 (M. d ) and <Pj be defined by © with scale factor 6j. 
Then for Algorithm^ there exist constants a 2 > 0, C5 > 0, Cq > such that 



I — 11 ^ n n+K(n- 

\U - U n ( K +l)\\L 2 (n) < ^4^1 



K 



Mliji(fl), 



K = 1,2, 



with ct\ given by (|20p . 



Proof Since there are K + 1 outer iterations (since the outer level index starts 
at 0) and n inner iterations, we have {K + l)n iterations in total, of which 
K (n — 1) + n iterations are with subsequently decreasing scale parameters 
for which we can use Theorem 15.21 The remaining K iterations involve the 
subsequent error estimation for K > and I = 1 since in this case, we have 
an increasing scale parameter for which Theorem 15.21 does not apply. With 
the proof of Theorem 15.21 we can derive a similar result with increasing scale 
parameters. In this case, we need to change the right hand side of ([H?]) and 
the following line to 



< 



and then we can define 



a 2 



1/2 



(28) 



Note that a 2 > u\ since fx < 1 by definition and the constants are all positive 
which will mean a lower rate of convergence as compared to Algorithm [1] The 
remainder of the proof follows the same steps as the proofs of Corollary 15.31 
and we leave the details to the reader. □ 

Note that the original justification for proposing this nested multiscale 
algorithm was that the errors from Algorithm [T] appeared to be dominated 
by a global behaviour, suggesting the need to go back and fit on a coarse set 
of centres with a large support. This is a similar idea used in the multigrid 
method in th e finite element l iterature (jBrenner fc Scottll2008l Chapter 6.3). 
As stated in (IFasshauerl 120071 Chapter 44.3), this additional outer iteration 
is known from Kaczmarz iteration, which i s frequently used in the multigrid 
literature as a smoother ( McCormick 19921 iKang fc Kwak 19971 ) ). 
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7 Numerical experiments 

In this section, we present the results from applying the multiscale and nested 
multiscalc algorithms to various PDEs. 



7.1 Multiscale Algorithm 

In this subsection we consider two PDEs, the first without Dirichlct boundary 
conditions and the second with Dirichlct boundary conditions. 

The first problem is the Helmholtz-like equation with natural boundary 
conditions: 

— Au + u — f in f2, 
d 

— — u = on df2. 
on 

We take ft = \— 1,1] 2 and f(x,y) = cos(7ra) cos(7rj/). The outer unit normal 
vector is denoted by n. The exact solution is given by 

cos(7ra;) cos(7rw) 
" (a ' g)= 2^ + 1 ■ 

We again use the C 6 Wendland radial basis function given by (fiM)) . Wc 
used five levels for the approximation, with equally spaced point sets at each 
level. The number of points, N, and the fill distances, h, are given in Tabled] 
We note that the fill distances decrease by almost exactly one half at each level 
and hence we select /i = ^. The L2 and errors and condition numbers (k) 



Level 


1 


2 


3 


4 


5 


N 


25 


81 


289 


1089 


4225 


h 


3.5e-l 


1.75c-l 


8.75e-2 


4.37e-2 


2.19e-2 



Table 1: The number of uniform points used at each level and the associated 
fill distance for the multiscale Galerkin approximation example 



of the stiffness matrix are given in Table [5J The Li error was estimated using 
Gaussian quadrature with a 300 x 300 tensor product grid of Gauss-Lobatto 
points and the L m error was estimated with the same tensor product grid. 
The second example uses the Poisson problem 

— Au = f in ft, 
u = on dfl. 

We take fl = [— 1, l] 2 and f(x,y) = sin(7ra;) cos (fy). The exact solution is 
given by 

sin(7ra;) cos (f y) 

u(x,y) = w . 

1.257H 
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Level 


1 


2 


3 


4 


5 


Sj 


2 


1 


0.5 


0.25 


0.125 


IMI2 


8.000e-4 


2.145e-4 


1.059e-4 


7.009e-5 


5.178e-5 


1 1 e j 1 1 00 


1.721e-3 


7.273e-4 


3.764e-4 


2.147c-4 


1.394e-4 




1.608c+3 


3.125e+3 


4.159c+3 


4.576e+3 


4.710e+3 



Table 2: The scaling factors, approximation errors and condition numbers of 
the stiffness matrices for the multiscale Galerkin algorithm with a Neumann 
boundary condition. 



We again use the C 6 Wendland function as the kernel, with the same 5 levels 
as for the first example. To verify that (fl"5"|) holds, we first check for the basis 
functions. For the boundary norm, since fi = [—1, 1] , without loss of gener- 
ality, we consider the x = — 1 boundary of the domain only. Then we have 
boundary integrals of the form ( Fa sshauer 20071 Appendix D) 



|V#i-n||i a(fln) 



—6 



dy 



22y (I6y 2 7y \ / y 
S 2 \ S 2 5 \ 5 



(by 



603969552384,5 



11305 



For the interior, we have 



|V<V 2 



L 2 (Q) 









2453<5 2 
4845 



r(^j(r) ) drd9 



22r /16r 2 
~P [~5 2 ~ 



d6, 



7r i 

T + 1 



r 



dr d9 



where O specifies the support of ^ in i7 and hence this last expression is finite 
and does not depend on S. 

In prac tical applications, we need t o select a value of (3n satisfying /3/v > 
2Cff/S. In lGriebel fc Schweitzer! ( 2002 ). it is proposed to estimate Cn/VS as 
the maximum eigenvalue of the generalised eigenvalue problem, 



Bv = ADv, 



where 



Bij = J (V<^ • n) (V# 3 - ■ n) , 

dO 



(29) 
(30) 
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and 

Dij = J (31) 

a 

where i and j run over the indices of all the radial basis functions with support 
overlapping the boundary. The extra calculation involved in this step is not 
significant since the entries of B are required for the construction of the stiff- 
ness matrix and the set of centres overlapping the boundary will generally be 
small compared to the entire set of centres. The maximum eigenvalue can also 
be efficie ntly computed with a simulta neous Ray leigh- quotient minimisation 
method (jLongsine fe McCormic3ll98Cll) . 



Level 


1 


2 


3 


4 


5 


Si 


2 


1 


0.5 


0.25 


0.125 


IMh 


8.125e-3 


1.451e-3 


3.229e-4 


8.217e-5 


2.219e-5 


1 1 e j 1 1 00 


1.057e-2 


2.350e-3 


6.504e-4 


1.963e-4 


5.962e-5 




5.633e+5 


l.OOlc+6 


8.056e+5 


4.572e+5 


2.374e+5 



Table 3: The scaling factors, approximation errors and condition numbers 
of the stiffness matrices for the multiscalc Galcrkin algorithm with Dirichlct 
boundary conditions 



The results a re in Table Bland s upport the theoretical findings above. We 
note that whilst IWendlandl (|l998d l did not find convergence after the third 
level with a similar algorithm, this may be due to the approximations used 
to calculate the integrals, rather than the algorithm itself. The potential for 
errors i n integration to affect the performa nce of Galerkin techniques are well 
known (|Strandl973LlBrenner fe Scottll2008l) . To estimate the integrals, we used 
the MATLAB functions quad2d and quad with an absolute tolerance value of 
le -10 . We also estimated the non-zero integration range both to speed up 
the calculations as well as to reduce numerical error which can result if for 
example, we integrate over the entire domain [—1, l] 2 whilst the function only 
has a very small support. 



7.2 Nested multiscalc algorithm 



In this subsection, we consider the same example as in Section 17.11 however 
now with Algorithm [2] with K = 2 and n = 2. We use the first two levels of the 
example described in Section 17.11 as the inner iteration. We also use the same 
kernel. Our choice of K and n leads to a 6 level algorithm. A similar example 
was considered in (| Wendlandl 1 1 998 a Section 5), however a lack of information 
regarding the exact approximation spaces used for the inner and outer level 
iterations means we have not been able to compare our results. The results 
from this 6 level nested algorithm are in Table HI 
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Level 


1 


2 


3 


4 


5 


6 


N 


25 


81 


25 


81 


25 


81 


IMI2 


8.000e-4 


2.145e-4 


2.045e-4 


2.088e-4 


1.991c-4 


2.034e-4 


\\ e j II 00 


1.721e-3 


7.274e-4 


4.325e-4 


7.002e-4 


4.165e-4 


6.742e-4 


Kj 


1.608e+3 


3.125e+3 


1.608e+3 


3.125e+3 


1.608c+3 


3.125e+3 



Table 4: The number of centres, scaling factors, approximation errors and 
condition numbers of the stiffness matrices for the nested multiscale Galerkin 
algorithm 



The results indicate erratic convergence and approximation errors far infe- 
rior to those using Algorithm [TJ This is not surprising since Theorem 18 . 2 1 indi- 
cates convergence of our approximation to the best approximation to u from 
V\ + V% whilst in Algorithm [JJ our approximation is formed from V\ + . . . + V5. 



8 Analysis of convergence 



In this section we will focus on esti mation of the con vergence and verify- 
ing approximation orders. Similarly to Wendland ( 2010t ). we will also rewrite 
the convergence results in terms of fill distances, which is t he usual form of 
convergence res ults for radial basis function approximations (IWendlandl2"(¥)5l 



FasshauerJl2Q07). 



8.1 Multiscalc Galerkin algorithm 

We consider Algorithm [TJ with hi = \i and /ij+i = [ihj. Since p. is a constant, 
we can rewrite (f2T)|) as 

ol\ = C 3 fl. 

Then with Corollarv l5.3l we have that 

l|e„||L 2 (r2) := II" - u«IU 2 (fi) < Chl- a \\u\\ H i {n) , (32) 

with 

a := -logc 3 /log^. (33) 

Hence we can either express our convergence in terms of an exponent of h n or 
equivalently a™. 

It is of interest that the error bounds do not depend on the kernel used for 
the approximation spaces. Typically with a kernel which generates H T (M. d ), 
we see error bounds proportional to h T . Since our kernel for the error analysis 
generates H 1 (M d ), we have b} n . Henceforth, we analyse the convergence in 
terms of o,\. We can calculate estimates of «i, which we denote by Si, as 
follows 

~ _ IKIUgog) 

ll e n-l||L 2 (r2) 
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C WF 


C b WF 


5l,2 


0.128 


0.268 


"1,3 


0.318 


0.494 


5l,4 


0.507 


0.662 


5l,5 


0.618 


0.739 



Table 5: The estimated convergence rates 5?i in using the results for the Li norm 
errors from the first example in Section 17.11 with the C 2 and C 6 Wendland 
functions. 

We can see that the estimated values of ot\ are higher with the C 6 Wend- 
land function, which indicates that we should not necessarily expect faster 
convergence with a smoother Wendland function and consequently we should 
not expect an error bound proportional to h T . 



8.2 Nested multiscale Galerkin algorithm 

In this subsection, we will focus on considering convergence of the nested 
multiscale Galerkin algorithm in terms of ot\ and oli- Note that a bound for 
the error at level n(K + 1) in terms of the fill distance at level n(K + 1), such 
as in (|52")l , will not be possible here because the fill distance at level n(K + 1) 
only depends on n and not on K. In other words, increasing K has no effect 
on the final fill distance. 

An additional benefit of considering the nested multiscale algorithm is more 
estimates of cti, particularly for repeated applications of the inner iterations 
(when K > 1). Table [5] gives the estimates of a.\ and ai from considering 
successive L2 norm error estimates in Section 17.21 Successive error estimates 
will be of the form 

ll e »llL 2 (J3) 

I|ei-i||i 2 (r2) 

By definition of our nested multiscale algorithm, for i = n + 1, In + 1, . . . we 
have an estimate for ai and in all other cases, an estimate for ol\. 

Tabic 1731 presents the estimated convergence rates OL\ >n and 52, n using the 
results for the L2 norm errors from the example in Section [7.21 



Level 




f*2,n 


2 


0.268 




3 




0.953 


4 


1.021 




5 




0.954 


6 


1.022 





Table 6: The estimated convergence rates ct\ in and «2,n using the results for 
the Z/2 norm errors from the example in Section 17.21 
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Interestingly, the difficulties with convergence are not due to ct2 which is 
seen to be less than 1 in all cases. This is empirical evidence of the effective- 
ness of the smoothing nature of the inner iterations, in that after the inner 
iterations, the errors are again of a global nature and hence a return to a 
coarse grid is justified. Convergence is affected by the repeated application of 
the inner iterations for which 5i,„ is always greater than 1. Empirically, this 
appears to suggest that the more localised features have already been captured 
in the approximate solution. 

This also implies that the angles between the approximation subspaces are 
close to zero since the linear convergence rate in Theorem 16. II can be bounded 
above by the angle between the subspaces. This is covered next. 



8.2.1 Angles between subspaces 



For our analysis of the convergence of the nested multiscale algorithm, we 
require the following definition of the angle between subspaces. 

Definition 8.1 Let "Hi and TA.2 be closed subspaces of a Hilbert space W with 
U := T-L\ (~1 %2- Then the angle a between Hi and Hi is given by 

cosa = sup {{u,v) : u € Hi n U ,v e^fl U and \\u\\, \\v\\ < l} . 

It is well known ( Smith et al. 1977 . Deutsch 199Clt) that Algorithm [5] con- 
verges linearly in the following sense. 

Theorem 8.2 Let u* be the best approximation to u from V\ + . . . + Vj with 
respect to \\ ■ |_ff 1 (fi)- Let u be the approximation from Algorithm^ Let a.j be 
the angle between Vj and (q—j +1 Vi. Then 

\\u* - u\\m(a) < c K \\u\\ H i {n) , 



where 



c 2 < 1 



II sin2 a r 

3=1 



This means that we need to estimate sin ctj to o btain upper bounds fo r the 
convergence rate. We follow a similar approach to Beatson et al. ( 200fj| ). We 
firstly define modified sets of centres as X\ = X\ and 



Xi = Xi\ |J X j; i > 2, 

3=1 

and the corresponding approximation spaces as 

Vi = span|<Pj.(-,x),x G X^ 
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Then we need to find the supremum of the inner product of u G Vi and 
v G Ai+i where 



A' 



A+i= |J Vj = spaii^* tfj (-,x i ),x i G (J X; 



J=i+1 



j=i+l 



with 



1. With the matrix K'f 12 ^ given by 



K 



Kff } = (# tf4 (-,x i ),# dj (.,x j )),x i eX i ,x i G |J X is 



and with coefficient vectors and v for u and u respectively, we seek the 
supremum of fiK^ 12 ^v. We also define matrices and K^ 2 ^ as 



K {1} 

v 



and 



A 



A 



K 



,{2} 



j(') x i)> ('i X j')) 5 X » ^ (J X 3^3 e U ^ 



j=i+l 



Let the Cholesky decomposition of be L^Li. This is well-defined since 

K' 1 ^ is strictly positive definite and symmetric. Then \\u\\ 2 = /i T L^Li/x and 
letting 7 1 = Li/i gives ||m|| 2 = 7^71 • We can follow a similar approach with 
K^ 2 } which gives \\v\\ 2 = 7^72 with = L^L^ and 7 2 = L2ZA However 

in this case since K^ 2 ^ is the union of radial basis functions with (possibly) 
different scaling factors, we cannot be sure that is positive definite. In 

our example, K^ 2 ^ was always positive definite and we do not dwell further on 
this. Sufficient conditions for an interpolation matrix cons tructed with several 
scaling factors to be positive definite can be found in (jBozzini et al. 2004. 
(11)) which al s o requ ires the Fourier transform of a Wcndland function from 
Chernih et al.l (2011) to com pute a lower bound on the minimum eigenvalue 
as given in (|Wendlandll2005l Theorem 12.3). 
Then we have 



(u,v) =fi T K^v 

= M T Lf (L^) T K {12} L 2 ~ 
= 7iM7 3 , 

with M := (Lj~ 1 ) T K^ 12 ^L 2 ~ 1 . The supremum of the inner product is given by 
the largest singular value of M. We denote this supremum by sincfj and the 
results with {T^}f =1 are in Table [7] We note that since X, C Xj, sina^ is a 
lower bound on sino^. We chose to estimate sin 5^ since by removing nested 
centres from later levels, we had less difficulties with singular matrices. 
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i 


1 


2 


3 


4 


sin cti 


9.849e-3 


2.684e-2 


4.153e-2 


6.987e-2 



Table 7: The estimates of sin on with the approximation spaces {Vi}f = 
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